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^^ ' The optical conductivity, o"(w), of the two dimensional one-band Hubbard model 
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is calculated at finite temperature using exact diagonalization techniques on finite 
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clusters. The in-plane d.c. resistivity, pab, is also evaluated. We find that at large 



'^ ' U/t and temperature T, pab is approximately linear with temperature, in reasonable 

^ ■ agreement with experiments on high-Tc superconductors. Moreover, we note that 

(t(w) displays charge excitations, a mid-infrared (MIR) band and a Drude peak, 

X' 

Vh ' also as observed experimentally. The combination of the Drude peak and the MIR 

oscillator strengths leads to a conductivity that decays slower than l/w^ at energies 
smaller than the insulator gap near half-filling. 
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Experimentally, it has been observed that the in-plane d.c. resistivity, pab-, of the hole- 
doped high temperature superconductors is linear with temperature when the hole doping 
fraction is optimal, i.e. when the critical temperature (Tc) is maximum. |1]] This sim- 
ple phenomenological law is still one of the most puzzling features of the normal state of 
the cuprates. A possible explanation of this behavior using the Bloch-Griineisen formula 
(which is based on electron-phonon scattering) seems unlikely, |^ and thus mechanisms 
based on scattering by spin fluctuations have been proposed. The a.c. conductivity, a{uj), 
also presents interesting features. A mid-infrared band (MIR) has been observed inside the 
charge-transfer gap of the insulating parent compound. In addition, at small frequency 
(relative to the gap) cr(c(j) decays as l/cj, instead of the more standard Drude behavior 
l/cu^. This effect can be phenomeno logically described by an energy-dependent lifetime 
t{(jj) ~ uj~^. Several theoretical mechanisms have been proposed to explain these features. 
1^ Most are based on simple mean-field solutions of electronic Hubbard or t-J like mod- 
els, but the validity of these approximate descriptions is unclear. An alternative approach 
involves direct numerical analysis of these models. Recently, there has been considerable 
progress in this approach and several studies of o-{uj) on finite clusters using exact diago- 
nalization techniques (at zero temperature) have been reported and compared to analytical 
approaches. PJ^ The presence of the MIR band has been explained as due to the consid- 
erable spectral weight located in the incoherent part of the hole spectral function, and the 
anomalous l/uo decay was attributed to a combination of oscillator strength between the 
MIR band and the zero frequency Drude peak at zero frequency in the metallic regime. Q 

In this paper we report a numerical study of the two dimensional (2D) one-band Hubbard 
model at finite temperature using the exact diagonalization approach on small clusters. 
Little work has been carried out previously at nonzero temperature using this technique, 
since the full set of eigenvalues and eigenvectors of the finite cluster is needed to determine 
thermal properties. This substantially increases the memory and CPU requirements relative 
to zero temperature properties. [0 Here we evaluate both a{ijj) and pab and compare the 
results with experiments; we find encouraging qualitative agreement. The calculation of 



transport properties of a weakly dissipative system in the context of many-body problems 
generally follows the Kubo formulation, which relates the conductivity to a current-current 
correlation function. This approach has been widely used in the context of Hubbard-like 
models to describe strongly correlated systems. The real part of the conductivity at finite 
temperature is given by 

a(a;) = 7rii^^^e-^^"|(n|j.|m)|25(u; + E„-Ej, (1) 

•^^ n,m 

where \n) is an eigenstate of the Hubbard Hamiltonian with eigenvalue En, Z is the partition 
function, j3 the inverse of the temperature, and j^ the current operator in the x-direction. 
The rest of the notation is standard, and details can be found in textbooks. |^ 

The diagonalization of the Hubbard model was carried out on small square clusters. In 
each subspace corresponding to a given set of quantum numbers (momentum, z-component of 
the total spin and parity under spin reversal) we computed all the eigenvalues and eigenvec- 
tors in two steps. First, the matrix was reduced to a tridiagonal form using the Householder 
algorithm. We then diagonalized the resulting matrix using a standard QL algorithm. Since 
in general we have to deal with complex hermitian matrices, we developed hermitian versions 
of the subroutines TRED2 and TQLI of the Numerical Recipes package. In principle, the 
total operation count for both subroutines scales as ~ N^, where Nh is the dimension of 
the matrix to be diagonalized. However, since the innermost loops could be vectorized, we 
found that the coefficient of the N^ was four orders of magnitude smaller than the coefficient 
of N^ term on a Cray YMP. The total CPU time required to diagonalize the largest ma- 
trix with Nh = 540 was approximately 23 seconds on a Cray YMP supercomputer and the 
total memory required for the diagonalization of a Nh x Nh matrix was 2 x Nh + 4 x Nh 
words. Both CPU-time and memory requirements compare well with similar subroutines 
in other packages such as IMSL or NAG. The calculation of cr(u;) itself was considerably 
more CPU time consuming than the diagonalization procedure. According to Eq. (1), the 
total operation count in a subspace of dimension Nh and for a fixed temperature scales as 
Nh. However, in the calculation of the matrix < n|jx|m > one can take advantage of its 



sparse nature, thus effectively reducing the dependence from N^ to N^. Moreover, a set of 
measurements at different temperatures could be done with almost the same CPU time as 
a single temperature by appropriately rearranging the loops and vectorizing the innermost 
one. 

Although the technique applied here works equally well for the Hubbard and t-J models, 
we have concentrated only on the Hubbard model, which possesses excitations across the gap 
that are important for comparison of cr(u;) with experiments. These calculations were carried 
out on small square clusters of eight and ten sites, similar to those used previously in the 
study of the Heisenberg model and other systems. 0,§] On finite systems it is important to 
choose the boundary conditions to appropriately minimize finite size effects. In the present 
study we decided to use antiperiodic boundary conditions (APBC). ||lT| In the 8-site cluster 
with APBC the non-interacting limit U/t = has 4 levels with energy -2t, and another 4 
levels with energy +2t. In the half-filled case, and with U/t > 0, two bands exist separated 
by a gap which grows as U/t increases. This behavior is also expected for the Hubbard model 
in the bulk limit. When holes are introduced, they are energetically favored to appear in the 
lower band. In contrast, the same cluster with periodic boundary conditions (PBC) contains 
6 levels of zero energy at U/t = 0, one state with energy -2t and another with energy +2t. 
The large zero-energy degeneracy appears to produce large finite-size artifacts in the PBC 
case at finite coupling, and for this reason APBC will be used in this paper. 

The application of Eq.(l) to Hubbard-like models involves some complications. One 
problem is that an isolated, finite system, such as the clusters analyzed in any computational 
study, cannot show resistive behavior. Thus, the resistivity of a metallic ground state at zero 
temperature must vanish, since a Drude-like weight D6{uj) appears at zero frequency in the 
conductivity. In fact, Eq.(l) implies that a 5-function is always present at zero frequency 
if any eigenstate of the Hamiltonian satisfies |(ra|jj.|r7i)| 7^ with En = Em- Therefore, 
electron-electron interactions (umklapp processes) are not sufficient to produce dissipation, 
and the d.c. resistivity of the Hubbard model is zero at all temperatures. [|T^ On finite 
systems with twisted BC the calculation of the Drude weight is carried out indirectly, using 



the two dimensional sum rule 

^°°dc..M = |^{-r), (2) 

where (T) is the thermal average of the kinetic energy operator. Assuming the existence of 
a contribution D6{uj) at zero energy, we obtain 

D (-f) 1 



where both terms on the r.h.s. of Eq.(3) can be calculated numerically. Following this proce- 
dure, the inverse of the Drude weight D~^ is plotted in Fig.l as a function of temperature for 
several couplings U/t, and at a filling of 6 electrons on the 8 site cluster {{n) = 0.75). D~^ is 
proportional to the d.c. resistivity, once a finite width is given to S{uj) to mimic dissipative 
processes not included in the Hamiltonian. It is interesting to note that for T > t and strong 
coupling, D^^ is approximately linear with temperature. This is in agreement with the pre- 
dictions of Rice and Zhang |13| for the large U/t limit. On reducing the Hubbard coupling 



U/t, we find that D^^ acquires curvature in T, and in the weak coupling region D~^ ~ T^. 
This behavior is consistent with the quadratic temperature dependence of the resistivity 
expected for a Fermi liquid. Now we will consider how our results can be compared with 
experiment, |^. First, note that for t = 0.4eV a temperature of It corresponds to approxi- 
mately 4600K which is much higher than the maximum experimental temperature for pab of 
~ 800K. In principle we should reduce the temperature in our cluster calculations for this 
comparison. Unfortunately, at experimentally relevant temperatures the finite size effects on 
the cluster are greatly increased, so that erratic behavior of the Drude weight as a function 
of (n) and U/t is observed. We estimate that for temperatures smaller than t/4 ^ 1200K 
our finite-cluster results are not representative of the bulk limit. It is thus more convenient 
to extrapolate the experimental results to higher temperatures since the slope dpab/dT is ac- 
curately known experimentally. These slopes (at the optimal doping concentration) are very 
similar among the different cuprates, and range from dpab/dT ^ IpQcm/K for Bi2Sr2Cu06 
(Bi2201) to dpab/dT ^ O.bpQcm/K for YBasCusOy-^ (Y123) with a T^ of 90K. The extrap- 
olated experimental results are shown in Fig.2 (dotted lines). The theoretical predictions 



obtained from the present cluster calculations are also shown in this figure (open squares 
and triangles), and were obtained by plotting D~^ times a parameter with units of fiQcm, 
which sets the relative scale between our calculations and experiment. (Physically this pa- 
rameter contains information about scattering processes not incorporated in the Hubbard 
model, so the overall normalization of our predictions is not determined and has been taken 
from experiment.) Given this freedom to fix the overall normalization, good agreement is 
observed between theoretical predictions and extrapolated experimental results for pab over 
the range of temperatures for which we consider the cluster results reliable. This encour- 
aging result suggests that the simple one-band Hubbard model can describe several normal 
state properties of the cuprates. There is a slight upward curvature in the results, which is 
not surprising since the experimentally measured pab is linear with temperature only at one 
particular density. [|1^ The filling fraction we have used, (n) = 0.75, may correspond to the 



slightly "overdoped" regime of the cuprates. flS 



Now let us analyze the a.c. conductivity. To allow a comparison with experiment we 
follow Imry []TB| and give each (5-function of Eq.(l) a finite width, to account for scattering 



through other processes not included in the model, such as phonons and disorder. This width 
e should be larger than the mean inter-level spacing in order to mimic a continuum of states, 
and thereby produce dissipation, e^^ can be considered to be a phenomenological relaxation 
time introduced to account for dissipative processes not included in the Hamiltonian. e is 
a free parameter in our study (in addition to the electronic density (n) and the coupling 
U/t of the Hubbard model) and we adopt e = 0.33 for the following discussion. In Fig.3, 
(y{uj) is shown for the 8 site cluster at several densities and couplings. The Drude peaks 
at zero frequency are incorporated in the plots. Fig. 3a shows the result for a temperature 
comparable to the antiferromagnetic exchange coupling J (T = 0.3125t), and for illustration 
we use U/t = 20 to enlarge the gap in the results. At half-filling (8 electrons), most of the 
spectral weight is located at cj > 5 (in t units), in other words these are charge excitations, 
as expected. As the density (n) is decreased, spectral weight is transferred from the high- 
frequency charge excitations to lower frequencies. A Drude peak is formed, and considerable 



weight appears within the insulating gap. (This may be associated with the MIR band 
observed in the cuprates as has been discussed extensively in the literature. |^) Fig. 3b shows 
the same cluster at T = 1.25t and coupling U/t = 8, which may be more representative of 
the cuprates. 0] Qualitatively, the behavior is similar to that found at lower temperatures 
and larger couplings, albeit with a smaller gap. Little sub-structure is observed in the 
spectrum. These results appear quite similar to the experimental observations of Uchida et 
al. on La2-xSrxCu04 (La214). Results for other high-Tc cuprates are very similar. Even the 
appearance of what Uchida et al P] called an "isosbestic" point (a point where conductivities 
for different densities cross) is reproduced in this figure. Finally, in the high temperature 
regime (T = 5t in Fig. 3c), the gap is completely filled at all densities, although a remnant 
of the upper Hubbard band can still be seen. The MIR band and Drude peak have merged 
into a single structure. In Fig. 3d, we show cr{uj) in the region 1 < tu < 5 for U/t = 8 
and T = 1.25t. For the case of < n >= 0.5, o"(u;) can be accurately described by a l/cu^ 
law, as expected for a conventional Fermi liquid. In contrast, for < n >= 0.75, cr{uj) has 
a much slower decay with u and can be fitted by the form (l/cj"), with a = 1.3 ± 0.2. 
Both forms are included in this figure for comparison. As can be seen in Fig. 3b for 7 
electrons (< n >= 0.875), cr{uj) has an even slower decay with u for u less than ^ 3. This 
anomalous frequency dependency of the conductivity has also been observed experimentally 
in La214. Moreover, a close correlation between the temperature dependence of pab and 
the frequency dependence of the scattering rate 1/t{uj) (seen in the frequency dependence 
of the conductivity) was observed; 1/t{uj) ~ u^'^ behavior in an overdoped sample of La214 
was recently reported, |T^ which is quite reminiscent of our power law fit. 



In summary, we have reported a numerical exact-diagonalization calculation of the op- 
tical conductivity of the two dimensional one-band Hubbard model. The d.c. resistivity 
shows linear behavior in T for T > t and large U/t, which presumably is also valid in the 
lower-temperature range J < T < t. For these temperatures our results compare well with 
experiment. The a.c. conductivity was also calculated and we have presented results at 
several temperatures. A MIR band is observed, together with a Drude peak and charge 



excitations. The combination of the Drude and MIR oscillator strengths leads to a conduc- 
tivity near half-filling that decays somewhat more slowly than l/a;^ at energies smaller than 
the insulator gap. 

Upon completion of this work, we learned of an independent study by Jaklic and 
Prelovsek (University of Ljubljana preprint, Dec. 1993) of the 2D t-J model at finite tem- 
perature using a new numerical method. Their results are qualitatively similar to our results 
for the Hubbard model. 
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FIGURES 
FIG. 1. The inverse of the Drude weight, D^^, obtained numerically on an 8-site cluster with 
APBC, as a function of temperature (in units of the hopping parameter t). The results are shown 
for couplings ranging from strong (U/t = 20) to weak coupling (U/t = 4). The filling fraction is 
shown in the figure. 

FIG. 2. The d.c. in-plane resistivity pab as a function of temperature. The dotted lines 
correspond to experimental results for Bi2201 and Y123, extrapolated to high temperatures com- 
parable to the hopping parameter (t is taken to be 0.4 eV (~ 4600K)). The squares and triangles 
are numerical results obtained on the 8-site cluster with APBC for U/t = 20 and filling fraction 
(n) = 0.75. D~^ was multiplied by a constant with units of inverse time to set the scale. The value 
of this constant was chosen independently for the two compounds. 

FIG. 3. (a) The real part of the optical conductivity as a function of frequency at different 
densities, for U/t = 20, e = 0.33, and T = 0.3125t; (b) as in (a) for U/t = 8 and T = 1.25t; (c) as 
in (b) for T = 5t; (d) (t{lo) at U/t = 8 and T = 1.25t in the interval 1 < w < 5, together with fits 
to l/oj^ (< n >= 0.5) and 1/lj^-^ (< n >= 0.25), indicated by dotted lines. 
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